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Penetration of Gamma Radiation From a Plane 
Monodirectional Oblique Source 1 

Martin J. Berger 

The penetration of gamma radiation from a plane monodirectional oblique source is 
solved by the moment method of Spencer and Fano. The novel feature of the treatment is 
the systematic investigation and exploitation of a function-fitting technique, which consists 
of approximating the desired radiation-flux function, f(x), by an expression of the form 
*2A n <j>(B n x), the parameters A n and B n being determined from the knowledge of the flux 
moments. 

A discussion is given of the criteria for the proper choice of the " weight function," 
<f>(x), and of the accuracy of the fitting procedure. The results of a sample calculation are 
presented for an 0.66-Mev source in water. 

1. Introduction 

One of the most effective techniques for calculating the diffusion of gamma raj^s in infinite 
media is the moment method [l], 2 i. e., the determination of the spatial and angular moments 
of the radiation flux and subsequent calculation of the flux-distribution function by moment- 
fitting. This technique has been successfully applied to the calculation of radiation from plane 
and point isotropic sources and from monodirectional plane sources emitting radiation per- 
pendicularly to the source plane [2]. The purpose of this work is to extend the method to the 
determination of the propagation of radiation from monodirectional oblique plane sources, i. e., 
sources emitting photons at a specified obliquity angle with respect to the source plane. In 
the problems treated previously, analytical asymptotic theory was available describing the 
radiation field at great distances from the source [3] ; this facilitated the reconstruction of the 
flux function from its moments. For the problem of the oblique source, reliable guidance in 
the form of a corresponding analytical asymptotic theory is lacking, so that increased care must 
be given to the moment-fitting procedures. 

Although specific application of the method described in this paper will be made to the 
diffusion of gamma radiation, the technique used should also be useful for determining the 
propagation of other types of radiation, e. g., neutrons, from plane oblique sources. 

We consider an infinite homogeneous scattering medium that contains, in the plane 2 = 0, a 
radiation source emitting 1 photon per square centimeter per second. The source photons are 
all emitted with energy E , and in directions inclined at an angle a with respect to the normal 
to the source plane. (The distribution of the azimuthal directions with which the photons 
are emitted may be considered arbitrary, being irrelevant for the propagation as a function 
of distance from the source plane.) Our task is to determine the energy flux, F (a) (E,z)dE, 
at a distance z from the source plane, due to photons in the energy range (E,E-\-dE) as a 
function of the obliquity a. The energy flux is defined as the flow of energy through an iso- 
tropic detector with unit cross-sectional area. The energy, E, will be expressed in million 
electron volts and the distance, z, in units of the mean free path of the primary radiation. 

Unscattered radiation will be treated separately because its direct calculation is trivial. 
In order to facilitate the moment-fitting procedures, separate consideration will also be given to 
singly-scattered radiation. This separation is desirable because its flux distribution differs 
considerably from that of the multiply-scattered radiation. 



1 This work was supported by the Office of Naval Research and the Armed Forces Special Weapons Project. 

2 Figures in brackets indicate the literature references at the end of this paper. 
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Section 2 discusses the moment-fitting procedures appropriate to the problem. In section 
3 specific application of these methods is made to radiation with a source energy of 0.66 Mev 
(cesium-137 radiation) in a scattering medium of water. Various energy spectra and space 
distributions are presented as functions of the obliquity angle. In section 4 the accuracy of 
the results is investigated by the comparison of various methods of moment-fitting and by 
comparisons with other related calculations and with an experiment. Details and subsidiary 
questions are given in appendixes. 

In appendix A, equations are derived for the calculation of the flux of singly-scattered 
radiation. In appendix B it is explained how the flux moments are calculated by the Spencer- 
Fano method; in particular, it will be pointed out how the moments required for the treatment 
of the monodirectional oblique plane source radiation can be obtained from a certain set of 
standard moments computed on the NBS automatic computer (SEAC). Appendix C presents 
details of the moment-fitting procedure. 

2. Moment-Fitting Methods 
2.1. Introductory Remarks 
We assume that the first eight spatial moments of the energy flux, 

z*(E,a)=f" n z n F^(E,z)dz, (1) 

(11=0,1, ... ,7), 

are known. (It is assumed that F U) stands for the flux of singly-scattered or multiply-scattered 
radiation, whichever we may want to determine. For details concerning the computation of 
the moments, see appendix B.) 

Although a well-behaved function is determined uniquely if all its moments are known 
the reconstruction of a function from only eight moments requires some finesse. It is desirable 
to incorporate prior knowledge concerning the flux that is available from related analytical 
calculations, physical intuition, or experience with similar problems, into a so-called weight 
function, which serves as the initial approximation. Successive improved approximations 
are then found by the use of the moments. 

In order to choose this weight function properly, we must give some thought to the ex- 
pected spatial dependence of the flux function. Assuming that 0<a<7r/2, i. e., that the 
photons are released by the source in the direction of positive z, we anticipate that the flux of 
scattered radiation will be a unimodal asymmetric function of z with a peak at a value of z 
somewhat larger than zero, and decreasing more rapidly on the negative than on the positive 
side of the peak. The asymmetry of the function will depend on the obliquity, being a maxi- 
mum when a=0, and zero when a=w/2. Furthermore, it can be expected that at large dis- 
tances from the source plane the flux has a dominantly exponential behavior as function of z. 

To find a weight function embodying these features is difficult. The task is made easier 
if the flux is split into a symmetric component f (a) (E, z) and an antisymmetric component 
g {a) (E,z), so that 

F^(E,z)=f^(E,z)+g^(E,z) (2a) 

F^(E,-z)=j^(E,z)-g^(E,z), (2b) 

where z denotes \z\. In this manner we can account for the expected peak without incorporat- 
ing this feature into jf (a) and g {a) , which may be chosen to be more or less smooth, monotonic 
asymptotically exponential functions. There is another advantage inherent in this splitting 
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technique. We can fit four even moments to determine f U) and four odd moments to deter- 
mine (f a) . This is much easier than fitting eight moments simultaneously. There is a slight 
disadvantage also; the flux cannot be determined very accurately far behind the source be- 
cause for large negative values of z the functions f U) and g (a) will have nearly the same value, 
and their difference will be rather inaccurate. This is not a serious drawback because the 
small amount of radiation far behind the source usually is not of much interest. 

2.2. Polynomial Expansions 

We shall now discuss the method of polynomial expansion used by Spencer and Fano 
[1] for the solution of the plane source with a=0 and explain why this approach cannot readily 
be carried over to the treatment of the source with a^O. 

Spencer and Fano also made use of the splitting technique and set 

p-WZ oo 

f m \E,z)=- ir T l c n (E)U n (wz) (3a) 

£ 71=0 

gW {Ej Z )J^— J2 d„(E) V n {wz), (3b) 

^ 71 = 

where z denotes \z\\ e~ wz /2 and wze~ wz /2 are the respective weight functions, 
and 

^=^(#,-0 v ^ (4b) 

(71=0,1,2, • • • ) 

are orthogonal polynomial systems associated with these weight functions. The coefficients 
c n (E) and d n (E) are related to the moments as follows, the nth coefficient involving only 
moments of order not exceeding n. 

C «(^=gQ)yW (5a) 

dm=%( n )yg^^W^ (5b) 

The validity of such a polynomial expansion depends on the correct choice of the exponent w 
in the weight functions. It can be shown by physical arguments and by related asymptotic 
penetration theory that a satisfactory choice for the plane source with obliquity zero is 
w ~ (up/no) , where ju is the linear narrow-beam absorption coefficient of the primary radiation 
and ix v that of the most penetrating component of the radiation. (This is usually, though not 
always, radiation with the primary energy.) With a correctly chosen parameter, the poly- 
nomial expansion will converge rapidly for wz up to —15, so that only a few terms requiring 
the use of the first few moments have to be fed into the expansion. 

The parameter w is difficult to determine for the plane oblique source. Clearly, the radia- 
tion component with the highest energy (usually the component with the smallest attenuation 
coefficient) will not govern the deep penetration because it must travel a distance approx- 
imately equal to z/cos a in order to get a distance z away from the source. On the other 
hand, radiation that has been turned parallel to the 2-axis travels away from the source plane 
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on the shortest possible path but has lost energy, so that its attenuation 'is increased. One 
would expect that the deep penetration is controlled by the small amount of the radiation 
that has been turned approximately parallel to the 2-axis in a number of collisions such that 
it has lost very little energy. It has been found by numerical experimentation that it is diffi- 
cult to obtain reliable polynomial expansions by adjusting the parameter w on the basis 
of convergence considerations. These difficulties can be overcome, however, by another 
approach [4]. 

2.3. Function Fitting 



If eight moments are available and an even-odd split is used, we set 

/ w (E, 2)=A 1 (E, a^zB^E, a)} +A 2 (E, a)<i>[zB 2 (E, a)} 
g^KE J z)=C 1 (E J amzD 1 (E J a)] + C 2 (E J amzD 2 (E f a)], 



(6a) 
(6b) 



where 4>(z) and \l/(z) are suitable weight functions with the desired asymptotic and symmetry 
properties : 

cp(z) = cp(-z) (7a) 



<p(z)=<p(z)e 
4*(z)=4iz)e 



*(g)=-*(-g). 

> for large 



(7b) 

(8a) 
(8b) 



where <p(z) and ip(z) are slowly varying functions of z. 

The quantities Aj and Cj are intensity (normalization) parameters, and the Bj and D } 
are scale parameters. The introduction of location parameters is not necessary because the 
peak in the flux function will be determined by the superposition of symmetric and antisym- 
metric solutions. The parameters Bj and D J7 analogous to the parameter w in the method of 
polynomial expansion, will be effectively determined from the moments themselves, so that no 
"inspired guess" is necessary. 

We require that 



HE } a)=T;Aj(E } a) P Z 2 ^[zBj(E,a)]dz 

j = l J -co 

T (#,a)=S^CE,<*) P z 2n ^[zDj(E,a)]dz. 

j=l J -oo 



(71=0,1,2,3) 



(9a) 



(9b) 



The solution of these equations is described in appendix C. The important aspect of the solu- 
tion is the fact that there are restrictions on the class of allowed weight functions <i>(z) and ^(2). 
In order to obtain the desired asymptotic behavior, it is required that the scale parameters 
B } and Dj are real-valued; otherwise, the weight functions would turn into oscillatory func- 
tions, which is inconsistent with the expected smooth behavior of the radiation density function. 
As shown in appendix C, two pairs of dimensionless parameters may be introduced to 
characterize the relation between the function to be fitted and the proposed weight functions: 



s=- 



m\ MoMA 



m m 4 



t = 



ml MIM, 



ra*m 6 



Ml 



mi 



(10) 



m\ M\M 7 



mim 7 



Ml 



y (ii) 



where M n =z n (E,a),m : 



J CO 



z 2n <p(z)dz, and m 2n+ i= 



■/; 



~2n+l 



\p(z)dz. 
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The weight function <p(z) belongs to the allowed class (i. e., the scale parameters Bj are real) 
if, and only if, the point (s,t) lies in an "allowed region" of the s-t plane, indicated by the shaded 
area in figure 1 and specified in appendix C. To determine whether ^(2) is allowable, we replace 
the point (s,t) by (s,t) and use the same diagram. 

It has been found in this investigation and in other related work, that considerable searching 
is often necessary to get an allowed weight function, but when there are several such functions 
they yield moment-fitting results that are in close agreement. Thus the function-fitting 
technique seems to operate in a sort of all-or-nothing fashion, and differs in this respect from 
the polynomial expansion method. 

The symmetric and asymmetric weight functions have been related as follows : 



rP(z) = Z<p(z). 



(12) 



This implies that at 2=0, the antisymmetric component of the scattered flux vanishes and the 
total scattered flux is continuous. 3 The choice of z, rather than some other antisymmetric 
function of z as the proportionality factor in eq (12) is somewhat arbitrary but recommends 
itself by its simplicity and is consistent with what is known about the analytical behavior ot the 
flux carried by photons with energies close to the source energy. 



3 The unscattered flux is obtained as space integral over the source function and the flux of the nth order of scattering as a space integral over 
the flux of the n — 1st order of scattering. The source function is a delta function 8(z) . The unscattered flux thus has a simple discontinuity at the 
origin, whereas the flux of higher orders of scattering is continuous. 




1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 



Figure 1. Diagram indicating whether a weight-function is "permissible" for function-fitting. 

The diagram should extend over the region 0<s< « , 0<t< oo , but only a portion of this region is shown. A weight-function 
is "permissible" if its (s, 0-point lies in the shaded region of the diagram. 
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3. Numerical Application: Radiation from a 0,66-Mev Plane-Oblique Source 

in Water 

A sample calculation, based on the methods outlined in section, 2 was carried out for a 
source of 0.66-Mev photons (radiation from cesium-137) in water. 

The flux moments needed as input data were calculated on the SEAC at energies ranging 
from 0.66 to 0.134 Mev. At lower energies, spectra were obtained through an extrapolation 
procedure based on two characteristic features of gamma-ray diffusion: 

(1) At very low energies an equilibrium spectrum prevails in the sense that the flux can 
(at least, approximately) be factored into the product of a space-dependent function and an 
energy-dependent function. This has been shown to be the case for radiation from point and 
plane isotropic sources [2, 3], and there is no reason to doubt that it will also hold for plane 
oblique sources. 

(2) The shape of the energy spectrum at low energies is nearly independent of the source 
energy. This is demonstrated by plotting in figure 2 the moment M (E) of the flux in water 
as a function of energy in the range from 0.134 Mev down to 0.020 Mev for source energies of 
0.255 Mev and 1.0 Mev. The two curves have been normalized in such a manner that they 
have the same ordinate at 0.134 Mev. It can be seen that they lie quite close together; the 
area under the corresponding curve for 0.66-Mev source energy, which must lie between them, 
can be estimated with an accuracy of 5 percent. (For the purpose of this estimate, only multi- 
ply-scattered radiation from the 0.255-Mev source should be considered because a 0.66-Mev 
source will produce only multiply-scattered radiation below 0.134 Mev.) 

Computations were carried out for the energies and angles of incidence indicated in table 1. 

Singly-scattered radiation was calculated directly according to the equations given in 
appendix A. Some moment calculations were also made and are discussed in section 4. 

Multiply-scattered radiation was treated by the moment-method according to the function- 
fitting procedure outlined in section 2. It is desirable that one can fit the flux for as many 
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Figure 2. Energy spectrum of scattered gamma radiation from 1.0- and 0.255-Mev sources in a 

scattering medium of water. 

These distributions were used to obtain by interpolation the low-energy part of the spectral distribution of radiation from a 
0.66-Mev source. 
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Table 1. Energies and obliquity angles for which the energy flux was calculated. 

[Asterisks indicate the values for which the multiply-scattered flux was determined by the method of function-fitting with the use of the weight 

function <p (z) = £s/2(z).] 



E 


a 


0° 


30° 


60° 


90° 


Mev 
0.660 
.503 
.325 
.259 
.215 
.184 
.150 
.134 


* 

* 


* 
* 
* 

* 


* 
* 


* 



energies and source obliquities as possible with the use of the same weight function. This not 
only inspires confidence that the choice of the weight function is correct, but more importantly, 
it makes possible smooth interpolation. A considerable number of weight functions were 
tried, and the results obtained with them are described in section 4. 
Preference was finally given to the function 



c P (z)=z 5/2 K 5/2 (z)=L 5/2 {z) 1 



(13) 



where K 5f2 is a Bessel function of the second kind with imaginary argument. This weight 
function was superior to other functions considered in that its region of applicability included 
most of the energies and obliquity angles of interest (see table 1). In the relatively rare 
instances of nonapplicability, interpolation with the aid of results obtained with some of the 
other weight functions was successfully used. In any case, it is shown in section 4.2 that the 
results obtained with different weight functions were generally in excellent agreement. This 
provides pragmatic justification for the moment-fitting technique. 4 

The results of the sample calculations for a 0.66-Mev plane monodirectional oblique source 
in water are illustrated in figures 3, 4, and 5. Figure 3 (a, b, c, d) shows energy spectra of 
the flux of scattered radiation at various distances from the source plane for obliquity angles 
a=0°, 30°, 60°, and 90°. 6 Figure 4 (a, b) shows the dependence of the total scattered flux 

/*0.66 

F ia) (E, z)dE on the source obliquity and the distance from the source plane. 

When the penetration of radiation as a function of the obliquity angle is known, one can, 
by superposition, determine the penetration of radiation from plane sources with arbitrary 
angular distribution. An example of such an application is the calculation of the ratio 



r(cos a', z)-- 



ri f0.66 

d cos a dEF^(E, 

Jcosa' */0 



z) 



Ji ro.m 

d cos a I dEF l «\E,z) 



(14) 



This is the fraction of the scattered integrated energy flux from an isotropic 0.66-Mev plane 
source that is due to source radiation emitted within a cone of opening a' around the normal 
to the source plane. This ratio is exhibited in figure 5 as a function of cos a' for various values 

of z. 



4 One suspects that there may be a deeper-lying reason for the wide applicability of Bessel functions of the type used in eq (13). According 
to Spencer (private communication) such functions have also proved themselves useful for fitting flux moments for a point monodirectional source. 

4 The flux is expressed as the flow of energy, in Mev per Mev of spectral energy per unit time, through a unit area of surface perpendicular 
to the direction of motion of the radiation, i. e., through an isotropic detector with an efficiency independent of the direction of the incident radi- 
ation. The normalization corresponds to a source emitting one photon with an energy of 0.66 Mev per unit time and unit area. Distances from 
the source plane are expressed in mean free paths (mfp) of the primary radiation. 
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Figure 3. Differential spectra F M (E, z) of the scattered gamma-ray flux ({Mev /cm 2 sec) I Mev) at various dis- 
tances z (mfp) from a plane monodirectional oblique source of 0.66- Mev radiation in a scattering medium of 
water. 

The ordinate scale corresponds to a source emitting 1 photon/cm 2 sec. 

a, Obliquity a=0°; b, obliquity a=30°; c, obliquity a=60°; d, obliquity a=90°. 
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Figure 4. Total scattered gamma-ray flux from a 0.66-Mev plane monodirectional oblique source in a scattering 

medium of water. 

The ordinate scale (in Mev/cm2 sec) corresponds to a source emitting 1 photon/cm2 sec. 
a, Flux versus obliquity; b, flux versus distance from source plane (mfp). 
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Figure 5. Fraction of the scattered energy flux from a 0.66-Mev plane iso- 
tropic source in water that is due to source radiation emitted within a cone of 
opening a' around the normal to the source plane. 
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4. Checks and Comparisons 
4.1. An Experimental Comparison 

Experimental evidence with which the theory could be compared is unfortunately quite 
scarce, in fact is nonexistent for the exact conditions of the present problem. The only relevant 
published work seems to be that of Kirn, Kennedy, and Wyckoff [5], who measured the pene- 
tration of 0.66-Mev gamma radiation from oblique plane sources through slabs of concrete. 

A comparison will be made by considering experimental and theoretical buildup factors. 
The energy buildup factor — a characteristic parameter widely used in shielding work — is 
defined as the ratio 

Total energy flux (scattered plus unscattered) 
Unscattered energy flux 

The fractional energy- transmission through a slab is the product of the buildup factor times 
the quantity exp (— z/cos a). 

The experiment of Kirn et al. and the calculation differ in a number of respects: (a) The 
calculation for an infinite medium allows for the possibility that radiation can penetrate beyond 
the detector and then be backscattered into it. In the transmission experiment, however, with 
the detector behind a concrete slab and with no further backing behind the detector, such 
backscattering is impossible. The buildup factor for a slab is therefore smaller than that for 
an infinite medium. In a Monte Carlo calculation [6] it has been found that the ratio of the 
two buildup factors becomes nearly constant for penetrations greater than 4 mfp. (b) It is 
only approximately true that concrete may be considered a water-equivalent medium (by 
expressing all distances in mfp of the primary radiation) . Actually, there is enough difference 
at low energies in the respective absorption coefficients (per electron) so that the buildup factors 
for concrete are appreciably lower than those for water. By referring to the extensive tabula- 
tions of Goldstein and Wilkins [2] it can be found that this reduction factor is a very slowly 
varying function of the penetration distance (except in the vicinity of the source). 

The reason for the approximate constancy of the buildup reductions due to effects (a) 
and (b) lies in the fact that both are reductions of the low-energy end of the gamma-ray spec- 
trum, but this part of the spectrum — as pointed out in section 3 — is an approximate equilibrium 
spectrum, i. e., rather independent of the depth of penetration. 

(c) The calculation assumes a monoenergetic source spectrum. In the experiment the 
source was physically large, and the source radiation was unfiltered, so that it presumably 
contained an appreciable amount of low-energy scattered radiation. Additional scattered 
radiation could have been introduced by the collimating arrangement. 

(d) The calculation is for the energy flux, whereas the experimental detector response is 
proportional to the energy flux for energies down to 200 kev only, but rises somewhat more 
steeply at lower energies. 

Figure 6 shows semilogarithmic plots of the experimental and theoretical buildup factors 
versus oblique penetration distance z/cos a (in mfp), for obliquity angles a=0° and a=60°. 
It can be seen that in both cases the experimental curves lie below the theoretical curves, and 
that for (z/cos a) gT3 the corresponding curves are almost parallel to each other, in conformity 
with the expectation of a constant buildup-factor reduction. For medium and deep penetration 
the theory therefore correctly predicts the relative attenuation of the energy flux at different 
depths. 6 



6 For a mere detailed quantitative comparison for obliquity a=0°, taking into account the boundary effect according to the Monte Carlo 
method, see [6]. 
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Figtjke 6. Theoretical and experimental buildup factors for radiation from a 0.66-Mev plane 
monodirectional oblique source, plotted versus oblique penetration distance. 

The theoretical curves pertain to an infinite scattering medium of water, the experimental curves to a plane-parallel 
concrete barrier. 



4.2. Comparison of Different Weight Functions and Moment-Fitting Procedures 

In the absence of much experimental evidence that would provide a check on the theory, 
one is forced to rely heavily on tests of internal consistency. Hence a good deal of numerical 
experimentation has been accomplished, employing various weight functions to be used in the 
method of function-fitting, including also polynomial expansions whenever possible. 

Listed below are the relevant properties of a few functions that were found to be useful as 
weight functions in fitting the flux distribution from monodirectional plane oblique sources. 
These functions were selected because they have the proper asymptotic behavior and are either 
available in tabulated form or easily evaluated directly. In the formulas below, z will be under- 
stood to mean \z\. 



Bessel Functions: 

<p(z)=L p (z) = z*K p (z), 
where K p is a Bessel function of the second kind with imaginary argument. 



™2n=r dzz 2n <p{z)=2 2n +vT Q+ra) r (!+tt+p)- 



(15) 



(16) 
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UmL p {z)=(^j n z*-^e->. (17) 



Up=q+i, (3=0,1,2, . . .), 

L^z)=($j^->±^m*^ (18) 

Generalized Exponential Integral: 

<p(z)=E p ( Z )=j i *-*'%■ (19) 

(For p=l, this is the ordinary exponential integral). 

Lim E p (z)=— (20) 

m - ( 4n > ! , ten 

Error Function: 

p(s)=l-erf(Vi)=l — 1= fV' 2 <ft, (22) 

V^r Jo 

Lim |~1 -erf (Vs)T=4=^ (23) 



Gaussian Weight Function: 



_ (4n+l)! 
M 2 „- 24B(2ri+1) , (24) 

p(s)=0(s)=^ (25) 

(2n)! ,„.. 

™ 2 »=24 (26) 

Table 2 gives a comparison of moment-fitting results obtained at an energy of 0.184 Mev 
for obliquity angles a=0° and 60°. Results for the symmetric and antisymmetric components 
are shown separately. For comparative purposes the normalization has been chosen so that 
the results obtained with the weight function L 5/2 are the standard of comparison and have the 
value unity. The over-all conclusion can be reached from an inspection of table 2 that the 
various methods are in good agreement; nonnegligible discrepancies occur only very close to 
or very far from the source. It is to be noted that at a=0° there is a good agreement with the 
well-established method of polynomial expansion into U- and F-poly normals. Even a Gauss- 
ian weight function yields essentially the same results as the others at distances not too far 
from the source. The results of table 2 are typical and have been confirmed by other similar 
calculations. 

It should be mentioned that function-fitting could also be carried out with the use of a 
larger number of moments than four by representing the flux by a functional series of more 
than two terms. If the use of the additional moments makes little difference, the fitting 
procedure is presumably accurate. Such a check would be equivalent to establishing the 
validity of a polynomial expansion by examining its convergence. In the sample calculation 
the number of available moments unfortunately was too small to do this. 
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Table 2. Comparison of the results of moment fitting with different weight functions 





Weight function (normal incidence: a=0°) 
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4.3. Other Comparisons 

The penetration of radiation from a plane isotropic source can be determined through the 
superposition of results for plane monodirectional oblique sources. Good agreement has been 
found between the results of such a superposition and a direct calculation based on a polynomial 
expansion with the use of the appropriate plane isotropic moments. The ratio of the respective 
results for the total scattered energy flux at various distances from the source is shown in table 3. 

As a further check, moment calculations were also carried out for singly-scattered radiation, 
which could be compared with the results of direct evaluations of the formulas of appendix A. 
Table 3 also contains sample results of such calculations; again, the agreement is good. It 
should be mentioned that the direct evaluation of the flux of singly-scattered radiation was 
rather laborious because of the necessary numerical integration. It would have been more 
expeditious, and would have involved only a minor loss of accuracy, to use moment-fitting 
making use of weight functions suggested by the exact analytical expressions or physical 
considerations. 

Table 3. Sample calculations to test the accuracy of the moment-fitting procedure 
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Ratio of results for singly-scattered flux from a monodirectional plane oblique source: 

A: Moment calculation/direct calculation («=0°, £=0.325 Mev; weight function Zs/2(z)). 
B: Moment calculation/direct calculation (a=90°; £=0.215 Mev; weight function Ei(z)). 

fir T0.66 

Ratio of results for scattered flux from plane isotropic source, I sin ada I deF^{E,z): 

C: Superposition of plane oblique source results/direct calculation based on plane isotropic moments. 
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5. Appendix A. Singly-Scattered Radiation 

Consider a photon emitted by a plane source at 2=0, with energy E Q and obliquity a. 
Assume that the photon travels to a depth s without interaction, undergoes a Compton scatter- 
ing between s and s-{-ds that changes its energy from E to a value in the range (E y E+dE), 
and crosses the 2-plane with obliquity co without further interaction. The probability for 
this contingency is 

PdEds== maE} exp j _M(go>_ M(£)( g -s) j s) cos 

cos a ( cos a cos co ) 

where ^(1?) is the narrow-beam linear total-absorption coefficient, K(E ,E) is the Klein-Nishina 

differential coefficient for Compton scattering with an energy change from E to E, and U(z) 

is the unit step function (Z7=l, if 2>0; £7=0, if s<0). According to the law of Compton 

scattering 

/ mc 2 mc 2 \ , . A mc 2 , mc 2 \ , AON 

cosco=cos a cos( 1 ^ — |- -^r- )+sm a sinl 1 =r-j—=- jcos <p, (A2) 

where mc 2 is the electron rest energy, and <p is an azimuthal angle distributed uniformly be- 
tween and 2ir. 

In order to obtain the energy flux F[ a) (E,z) of singly-scattered radiation, one must multiply 
eq (Al) by E, integrate over all values of s between and o° (assuming that 0<a<7r/2), 
multiply the result by sec co 7 , and average over all possible azimuthal angles <p. With distances 
expressed in mean free paths of the primary radiation, the final result for the singly-scattered 
flux can be written 



where 



F["\E, z)= E ff >E) e _ u/0Ma) 1 C ir d ^ H( E a) (A3) 

fjL{rLo) cos a atJ 

H(2,E,a)=^^ {(e zt -l)U(cos co)-C7(-cos co)} (z>0) 

e zt U(~cos co) (g<0). (A4) 



t=- 



t COS CO 

1 m(JE) 



cos a /jl(Eo) cos co 
The integral in eq (A3) has been evaluated numerically, except when a=0°; in this case 

(TflC Trie \ 
1 — ^"+£ r ~r ^ 5 ) 

and the integration over <p is trivial. 

6. Appendix B. Calculation of the Spatial Flux Moments 

The following brief discussion of the moment equations is presented to make this paper 
self-contained. It will also bring out some minor points in which the procedure followed here 
differs in details from that given in the basic reference [1]. These are the use of an energy flux 
instead of a number flux, the separation of unscattered and singly-scattered flux, and the 
extension to plane sources with arbitrary angular distributions with cylindrical symmetry. 



7 "Without the factor sec a one would obtain an energy current across the z-plane instead of the energy flux through surface perpendicular to 
the direction of propagation of the radiation. 
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Let F k (z, cos 6, E) be the energy flux of photons that have been scattered k times. In terms 
of the Ffc's, the transport equation has the form 



-^ cos d+vF =S(z,E, cos 6) 



^fcosd+»F k =j E dE'~K{E',E). 

dip' d cos 0' 5- 5 ( 1 — p-+-pr— cos cos 0'- 
sin sin 0' cos <p' J F k _ x (£>0) 

where 5 is the Dirac delta function, and S is the source function. 
We shall consider a set of elementary source functions 



(Bl) 



S L =E S(E-E )P L (cos 0)5(2) 



2L+1 
4tt 



£=0,1,2,..., 



(B2) 



where P L is a Legendre polynomial. These elementary sources have no direct physical 
meaning except when Z=0 but are convenient mathematical fictions. 

We shall indicate the source type from which a given flux is derived by an index L. 

We expand F kL into a Legendre series in the angular variable and take spatial moments of 
the Legendre expansion coefficients; i. e., we introduce the transform 



G kL (n,l,E)= I dz\ d cos 6z n Pi (cos 0) F kL (z, cos 6,E) 

J-oo J-l 



and a similar transform for the source function 

T L (n,l,E) 



■/>£• 



d cos dz n Pi (cos 6)S L (z, cos 0,2?)=— E 5(E—E ) d n0 d L 

JiTT 



(B3) 



(B4) 



It follows from the transport eq (Bl) and sundry properties of Legendre polynomials, that the 
transforms G kL obey the equation 



(A+j0<%l=2* 

(A+»)G kL =fyE^-, K(E',E)P t (l -^+|f) &_ M 

where A is an operator defined by the relation 



(k=0) 
(*>0), 



AflS B (n,W=ngtl.flJB(n-l > iH-l^)+n 2 q^ ftt(n-l,Z-l^). 



(B5) 



(B6) 



The transform of the multiply scattered flux is found by summation 

H^E^itG^nlE). (B7) 

k=2 

It follows from eq (B5) and (B7) that 

(A+ l ,)H L (n,l > E)=jyE'§ 7 K(E\E)P l (^-^+^Hdn^E')+Tt(n,l > E), (B8) 
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where the source function 

Tt{n,l,E)=£°dE'§ 1 k{E',E)P l (i-^+^f) G 1L {n,l,E'). (B9) 

The calculation thus proceeds as follows: One computes G QL and G 1L with the use of eq (B5), 
then the source function according to eq (B9), and finally H L according to eq (B7). The equa- 
tion for Gql is algebraic; all others are linear integral equations of the Vol terra type, and are 
solved numerically on the NBS automatic computer (SEAC). The equations are interlinked, 
but the nature of the coupling operator A is such that they need not be solved simultaneously. 
Indicated below are the chains of equations that must be solved to calculate H L (n,l,E) or 
G kL (n,l } E) ; for the sake of brevity only the (n,l) values are shown. 
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For example, in order to compute Hq(4, 0, E), one must first calculate in turn the quantities 
#o(0, 0, E), flo(l, 1, E), H (2, 2, E), and # (3, 1, E). 
It can be shown that 

H L (n,l,E)=G kL (n,l,E)=0 (B10) 

unless n—\L—l\ is a nonnegative even integer. 

Now consider a monoenergetic plane source with arbitrary distribution, expanded into a 
Legendre series: 

oo 07"J_1 

S(z, cos d,E)=E 5(E-E )5(z)^^p^ AJ> L (cos 0). (Bll) 

L = 47T 

It follows from the linearity of the transport equation that the moments of the flux for such a 
source can be obtained through the addition of moments for "elementary" sources. Taking 
into account (B10), one finds that 
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z>"(E)=T,A 2L Hu(2n,0,E) 



i=0 



z*»+'(E)=J2A 2L+1 H 2L+1 (2n+l,0,E). 



(B12) 



It is remarkable that in order to determine the first n even (odd) moments, one need specify 
the source only by its first n even (odd) expansion coefficients A L . 

Finally, we note that for a monodirectional source with angular part 8 (cos d— cos a), 

A L =P L (cos a). (B13) 

7. Appendix C. Details of the Function-Fitting Procedure 

For the solution of the moment eqs (9a, b) it is convenient to adopt the following 
notation: 

Mo 



A.j=oljBj 

R rm„A 



vMi 



3 \jn % M l '} j 



(C2) 



0=1,2) 

Equations (9a, b) can then be written in the form 



= 1,2) 



a 1 /3i + a 2 |S 2 =l 



(C3) 



Ti+T 2 = l 
7,§ 1 +T 2 5 2 =1 
7i5f+r 2 5l=S 

y 1 s\+y 2 dl=t 



Y (C4) 



Equations (C3) and (C4) being formally the same, further discussion is confined to (C3). 
A solution can be obtained by a method given by Chandrasekhar [7]. Let 



t-s 2 



s-t 



Wo= 



^ (C5) .,=;— (C6) 

It follows from (C3) that fr and /3 2 are the roots of the quadratic equation 

|8 2 +«ii8+wo=0, 



(C7) 



and that 



—fcs (C8) 



«i=l — <x 2 (C9) 



that 



Bi and B 2 are real-valued if and only if p l and fi 2 are real and positive. This requires 



coi<0 

co >0 



i 



(CIO) 



u>\— 4co >0 > / 
Equations (C5), (C6), and (CIO) imply the following conditions on 5 and t: 
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IfO<s<l,0<£<l, the allowed region lies between the parabola s 2 =t, and the curve whose 
parametric representation is 

S= I {S+Qu~u 2 -[(3+Qu-u 2 ) 2 ~^u] 1/2 } 1 

s l y (en) 

t=us J 

(0<u<l) 

If s> 1 and t> 1, the allowed region lies between the line 8=1 and the parabola s 2 =t. 



The author expresses his gratitude to L. V. Spencer and U. Fano for helpful discussions and 
suggestions. He also is indebted to J. Doggett and J. Hubbell for their help with the numerical 
calculations. 
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